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Abstract 

We propose an approximation to general relativity that captures the main 
gravitational effects of dynamical importance in supernovae. The conceptual 
link between this formalism and the Newtonian limit is such that it could 
likely be implemented in existing multidimensional Newtonian gravitational 
hydrodynamics codes employing a Poisson solver. As a test of the formalism's 
utility, we display results for rapidly rotating (and therefore highly deformed) 
neutron stars. 
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I. INTRODUCTION 



The collapsed cores of massive stars are relativistic bodies. In addition to relativistic 
effects, convection and rapid rotation may play important roles in core collapse supernovae, 
making their accurate simulation a three-dimensional (3D) endeavor. Because accurate and 
stable numerical solutions to 3D relativistic problems are difficult to achieve, an approxi- 
mation that captures the relativistic phenomena most relevant to supernovae is desirable. 
(More detailed discussion of these points, together with references, are given in the following 
paragraphs.) 

While relativistic treatments would be necessary to follow the detailed processes leading 
to black hole formation, even collapse events leading to neutron stars ought to be treated 
relativistically, particularly at late times. Denoting the gravitational potential by $, the 
Newtonian limit is obtained from the Einstein equations using the metric 

rfs^ = -(l + 2$)rft^ + (l-2$)rfx2, (1) 

provided $ <^ 1 and velocities are much less than the speed of light. Taking M and R to 
be the mass and radius of the collapsed core, the gravitational potential is expected to be 
of order 

R \IA MqJ VIO km/ ' ^ ^ 

while infall and/or ejecta velocities of order 

might be encountered. If some collapsed cores are born with rotation periods comparable 
to observed millisecond pulsars, equatorial velocities would be of order 

„_.„ = M = 0.2l(I-)"'(^). (4) 

Crudely taking the core as a cold degenerate gas of nucleons of mass tub and Fermi momen- 
tum pfi nucleon velocities of order 

clearly indicate that pressure and internal energy density cannot be neglected as gravita- 
tional sources. The numerical coefficients in equations (0-|^) may seem to overstate the case 
for relativity immediately after core bounce, when a less massive inner core has a larger 
radius. They are perhaps more relevant numbers for late times, when the final explosion 
energy and remnant nature (neutron star vs. black hole) are determined. However, even 
at earlier times relativity cannot be ignored: detailed comparisons of Newtonian and rel- 
ativistic simulations in spherical symmetry show more compact cores and higher neutrino 
luminosities and average energies in relativistic treatments 0,^. Because of indications that 
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small (several percent) variations in, for example, neutrino luminosities and shock stagna- 
tion radius can make the difference between a successful (exploding) and a failed supernova 

it is clear that even modest relativistic effects comprise an indispensible component of 
reahsm in supernova studies. 

The multidimensional nature of supernovae must also be recognized as an important 
aspect of realism. Various observations, especially data from SN 1987A, point to the as- 
phericity of supernova explosions (see e.g. for an overview, and for a recent polarimetry 
analysis of several supernovae). Convection — either deep in the core and driven by a lepton 
gradient [^] or doubly-diffusive phenomena ("neutron fingers") [[7|, or in the more tenuous 
region between the neutrinosphere and the stalled shock, driven by an entropy gradient 

— has for some time been suggested as a means of boosting neutrino luminosities or 
increasing the efficiency of heating the material behind the stalled shock. The results of 
various simulations in up to two dimensions (2D) differ as to the significance of convective 
effects [fO-|T^,IT],HOJT3] ; clearly further study is necessary. Beyond convection, a pioneering 



2D simulation studied jet production by magnetic fields in a supernova context [jT5[, and 
simulations of explosions from collapsed stars with phenomenologically introduced jets have 
also been performed in 2D [p!^-[T8|. 

These 2D simulations have yielded interesting results, but ultimately consideration of the 
third dimension will be necessary. Based on an initial exploration of 3D effects, it has been 
reported @] that the sizes of convective cells in 3D simulations are about half as large as in 
2D simulations. Moreover, rapid rotation can significantly affect the strength and spatial 
distribution of convection [|19|. Detailed studies of magnetic field generation, jet formation, 
and neutron star kicks also invite 3D treatments. 

Including all the physics necessary for realism in a supernova simulation is a daunting 
task. The multidimensional simulations mentioned above, which had simplified neutrino 
transport, taxed the computational resources of their time; the same is true of recent sim- 
ulations involving Boltzmann transport in spherical symmetry [^,^,^. Adding general 
relativity to the list of desired physics makes things all the more challenging. While nu- 
merical relativity has been successful in spherical and axisymmetric cases, "...in the general 
three-dimensional (3D) case which is needed for the simulation of realistic astrophysical sys- 
tems, it has not been possible to obtain stable and accurate evolutions...," and it is argued 
that the difficulties are more fundamental than insufficient resolution [E2l. While the diffi- 



culties with 3D numerical relativity are beginning to be overcome [^,^, the computations 
are "prohibitively slow" |^ when the central gravitating body has a modest compactness 
{M/R ^ 0.15). Hence it will be some time before supernova simulations with sophisti- 
cated microphysics and transport can cover the full process of collapse, explosion, and wind 
formation over hundreds of milliseconds in a fully relativistic context. 

In order to overcome the difficulties associated with 3D general relativistic simulations, 
and to save resources for 3D hydrodynamics and accurate neutrino transport, an approxi- 
mate treatment of gravity that captures the phenomena of dynamical importance in super- 
novae would be desirable. The list of new gravitational phenomena introduced by relativity 
includes the nonlinearity of the gravitational field, the inclusion of all forms of energy and 
stress as sources, and gravitational waves. The first two of these effects are of dynamical 
importance in supernovae, while gravitational waves will probably not exert a strong back- 
reaction (unless bar or breakup instabilities of some sort become operative in the core). 
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Given the effects we wish to capture, what sort of approximation might we try? A common 
procedure in the hterature has been to assume a spherical mass distribution in implementing 
gravity, which allows a relativistic treatment. However, for rapidly rotating progenitors this 
introduces errors of about 10% during collapse and bounce, with errors rising at later times 
as the nascent neutron star cools |T^. A multidimensional approach that captures rela- 
tivistic phenomena could reduce these inaccuracies. A 3D post-Newtonian approach might 
be considered, but the post-Newtonian limit is known to be inadequate for quantitative 
determinations of neutron star properties. The case of spherical neutron stars constructed 
with a polytropic equation of state (EOS) with adiabatic index F = 2 (a common proxy for 
more realistic dense matter EOSs) is illustrative: when the polytropic constant is chosen to 
make the maximum mass in a fully relativistic calculation equal to ~ 1.6 M©, the maximum 
mass determined by first-order post-Newtonian calculations is ~ 2.5 Mq — an error of close 



to 60% 



We conclude that an approach that both probes the nonlinear nature of gravity more 
deeply, and does so in multidimensions, would be desirable in the supernova context. A 
method that could be incorporated into existing Newtonian hydrodynamics codes would be 
even more useful. We describe such an approximation — which we call "NewtonPlus" — in 
§2, displaying a full set of Einstein equations in order to see where inconsistencies arise and 
how serious they are in the supernova environment. Hydrodynamics equations for the radial 
direction are presented in §3 in order to argue that the approach could be implemented in 
existing codes. In §4 we display results for rapidly rotating (and therefore highly deformed) 
stars, which show that this simple "NewtonPlus" approach to gravity is indeed a significant 
improvement over the Newtonian limit. Concluding remarks, including comments about 
certain approximate approaches to gravity employed in binary neutron star calculations, 
are given in §5. An appendix shows how the static limit of our formulation relates to the 
familiar equations of stellar structure in Schwarzschild coordinates. 



II. EINSTEIN EQUATIONS IN THE NEWTONPLUS APPROXIMATION 

As mentioned in §1, use of the metric of Eq. ([1|) in the Einstein equations yields the 
Newtonian limit, provided the gravitational potential $ -C 1 and velocities (including mi- 
croscopic velocities, large values of which lead to significant stresses) are much less than the 
speed of light. In order to capture the nonlinearity of gravity, the significance of stresses 
as gravitational sources, and relativistic fluid velocities, we propose the use of the following 
metric: 

ds^ = -e^'^+^'dt' + e-2*(t/r2 + r^d^""), (6) 

where dVl? = dO"^ + r^sin^^c?^^. In comparison with Eq. (0), the "linearized" metric 
functions have been promoted to full exponentials, and a second metric function, 5, has 
been added. Eq. (§) will then reduce to the Newtonian case if 5 — 0; we shall see that this 
is in fact the case under conditions in which the Newtonian limit is valid. This provides 
a tight conceptual link with the Newtonian limit. Since it has two independent metric 
functions, this "NewtonPlus" metric also provides an exact solution in spherical symmetry. 
(For the static case, the Appendix shows the connection between our formulation and the 
familiar equations in Schwarzschild coordinates.) 
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A convenient formulation of the Einstein equations is the (3+1) formahsm, in which 
spacetime is fohated into spacehke shces labeled by a time coordinate ||27H30|. The metric 
can be expressed 



ds'^ = _ p.pi)dt'^ + 2pi dt dx' + 7^- dx'dx^ , (7) 

where the latin indices run over the three spatial coordinates. As originally conceptualized, 
the four quantities a (the lapse function) and Pi (the shift vector) could be chosen at will as a 
"gauge choice" or as "coordinate conditions," corresponding to the invariance of general rel- 
ativity under coordinate transformations. Four constraint equations — the so-called Hamil- 
tonian constraint and momentum constraints — would be solved to provide self-consistent 
data on the initial spacelike slice. The six degrees of freedom to evolve in time would then 
be the jij, determined by equations second order in space and time. Alternatively, one could 
evolve equations first order in time for 'jij and Kij; the latter is the extrinsic curvature ten- 
sor, which describes the embedding of the spacelike slices in spacetime. As an alternative to 
prescribing a and f3i, one can, for example, place conditions on Kij; then a and /3j become 
quantities for which solutions must be found. 

It is apparent that the NewtonPlus metric of Eq. (^ contains only two of the six degrees 
of freedom that should be present. This means that examination of a complete set of Einstein 
equations should reveal inconsistencies; we here explore these and discuss their seriousness 
in the supernova environment. 

We begin with the Hamiltonian constraint. This yields 

= Ane-^^E + -{d<l>f - -e-^''-^\dt<!>f . (8) 

In this expression is the usual 3D flat-space Laplacian. The energy density as viewed 
by an "Eulerian" observer (i.e., one whose 4- velocity is orthogonal to the spacelike slices, 
having covariant components = {—a, 0, 0, 0), where a is the lapse function) is denoted by 
E = T^'^nfj,n^, where T'^'^ is the stress-energy tensor. For a perfect fluid, E = T'^{p + p) — p, 
where F = (1 — 11^)^^/^, f is the magnitude of the physical fluid velocity as measured by an 
Eulerian observer, and p and p are respectively the total energy density and pressure in the 
fluid rest frame. We have employed the notation 

dX dY = drX drY + \ deX deY + , ^ d^X d^Y. (9) 

As expected from the conceptual link between the Newtonian and NewtonPlus metrics, Eq. 
(§) identifies $ as a glorified gravitational potential. In addition to the rest energy, the source 
includes internal and kinetic energies and pressure, all boosted by nonlinear contributions 
from $ itself. 

Next we turn to the momentum constraints. These can be summarized by 

V(e-*"^9t<l') = -47re-*s. (10) 

The three momentum constraint equations are obtained by reading this as a vector equation 
in orthonormal spherical coordinates. For a perfect fluid, s = F^(p + p)v, where v is the 
physical velocity measured by an Eulerian observer. Equation (0) is the first challenge to the 
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consistency of the use of the NewtonPlus metric: Since the curl of any gradient vanishes, 
equation (10) will only have a solution if the fluid flow is such that V x (e~*s) = 0. 
This condition is satisfied in spherical symmetry, but it typically will not be true when 
convection is present. For the use of the NewtonPlus metric to be meaningful, there are two 
possibilities. One is that the transverse (or solenoidal) portion of e~*s is small compared 
with the longitudinal (or irrotational) part, that is, V x (e~*s) ^ V- (e~*s); this will obtain 
if, for example, radial flows dominate convective effects. In that case, e~*s is determined by 
the Poisson equation 

V\e~'''^dt^) = -47rV ■ (e-*s). (11) 

The second possibility — the one we expect to follow in practice — is to neglect 5^$. This is 
reasonable in supernova environment: In contrast to the binary merger problem, one does 
not have entire stars moving at relativistic velocities. In the absence of bar or breakup in- 
stabilities, there is a quasispherical core, with asphericities due to rotation and convection. 
Aside from secular changes due to viscosity, rotation alone can be considered relatively sta- 
tionary, not contributing greatly to explicit time variation of $. Core convection may occur, 
but probably will involve nonrelativistic velocities. There may be relativistic radial infall 
velocities, as seen at late times in failed explosions |0, and also relativistic outflow velocities 
in jets but these situations involve matter at low density in comparison with the core. 



Hence in the various phenomena occuring in the supernova environment there is reason to 
suppose that either the densities or velocities are such that s = r^(p + p)v is everywhere 
small enough for dt^ to be neglected altogether. Heuristically, ignoring explicit time deriva- 
tives in determining $ based on this physical reasoning is consistent with the expectation 
that the back- reaction due to gravitational waves will not be dynamically important. 
Next, we consider the evolution equations of 

dtjij = -2aK,, + 7,fcA/5' + %kD,l3\ (12) 

where D represents a covariant derivative with respect to the 3- metric 7jj. For the Newton- 
Plus metric, the extrinsic curvature turns out to be 

=e-*-'(9i$)7.,. (13) 

Referring back to Eqs. (^) and (0), it is easy to see that the vanishing shift vector conspires 
with the explicit form of the extrinsic curvature tensor to make Eq. (|T^) an identity leading 
to no new information. (This is the case even if we choose not to neglect dt^.) 

Finally, we come to the evolution equations for Ky. In considering these equations we 
will take (9($ = in accordance with the physical reasoning discussed in connection with the 
momentum constraints. From Eq. ([T3|), we see that our neglect of explicit time derivatives 
of $ implies a vanishing extrinsic curvature tensor. Nevertheless, the evolution equations for 
Kij provide nontrivial conditions, and will lead us to an equation for 6. Various combinations 
of the Einstein equations could be employed to give an equation for 6, and we have explored 
some of these in numerical computations of stationary rotating stars. Because of numerical 
convergence issues discussed at the end of this section, our experience shows that convergence 
is achieved when the equation for 6 is derived from a combination of the individual evolution 
equations for Kij] the "lapse equation," which is the trace of the evolution equations for 
Kij] and the Hamiltonian constraint. 
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The equations for S obtained in this manner are as follows. Combining the evolution 
equation of K^^ with the lapse equation and the Hamiltonian constraint yields 

^ ^d'P^f - T^-AdeSf, (14) 



2r2sin[^]^' ' 2r2 

while combining the K^g evolution equation with the lapse equation and Hamiltonian con- 
straint yields 



r 



r"^tant^ r'^sm 

r^sm 



A$9^5-(9,5)^-^-^(9<^5)^ (15) 



r^sm r^sm 6 

and similar use of the K"^^ evolution equation gives 

In these equations, S^y = T^^h^ph^^, where the spacelike projection tensor is defined by 
hp,u = dfiu + nf^n^. For a perfect fluid, S^^ = {E + p){v^Y + p, where is the component 
of the fluid's physical velocity in the 6 direction as measured by an Eulerian observer, 
and S'^^ is given by a similar expression with replaced by v'^. Given the fact that we 
expect inconsistencies in the Einstein equations due to the reduced number of degrees of 
freedom that we keep, it seems likely that Eqs. (|T^-(|TB|) for 6 are inconsistent, though it 
is not obvious (to us) how to prove this rigorously. However, we note that stresses (e.g., 
pressure) constitute the primary source terms in the equations determining 6, confirming the 
expectation expressed previously, that S should vanish as the Newtonian limit is approached. 
The observation that 5 will only be appreciable at the highest densities, where pressure begins 
to make a nontrivial contribution in comparison with energy density, suggests a reasonable 
path forward. In typical cases it is expected that the deepest portion of the core will be 
roughly spherical, even if rapid rotation causes an equatorial bulge of lower density material 
(e.g.. Fig. 2 of Ref. |^). If this is the case, neglecting the angular derivatives of S is 
justified, removing many of the apparent inconsistencies in the three equations for 6. The 
remaining discrepancies have to do with angular derivatives of $ and particular components 
of the stress tensor that appear in Eqs. (P!^)-([T^). As it happens, if one adds (or averages) 
Eqs. (|l^)-(0), these remaining discrepancies disappear. Hence the equation we shall use to 
determine 6 is 

drdr6 + ^dr6 = 47re-2*(5% + S'^^) - [5,($ + 5)f. (17) 
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The source on the right-hand side is to be angle-averaged in solving for S. Recalling that 
S^g = S'^^ in spherical symmetry, Eq. ([TtD is the equation for 5 obtained in the spherical 
case. 

It should be straightforward to include the solution of $ and, if desired, 5, in existing 
multidimensional gravitational hydrodynamics codes. The solution of $ would make use of 
the Poisson solver normally used to solve for the Newtonian gravitational potential, the only 
difference being that one would have to iterate on Eq. (with the dt^ term dropped) to 
get a self-consistent The simplest approximation would be to simply solve for $ in this 
manner, and ignore 5 altogether. The next level of approximation would involve solving 
Eq. (p!?!) for 5, but ignoring 5 on the right hand side. Since 5 is already something of a 
correction, 5 appearing on the right hand side is essentially a "correction to the correction." 
If desired, however, Eq. ( p!7| ) could be solved as it stands, with iteration being required. In 
any case, solving Eq. (0), while requiring only a couple of angular and radial integrations 
over the computational domain, is a bit subtle. This is because the Green's function of the 
operator on the left-hand side is a logarithm, which must be delicately handled in order to 
get a vanishing boundary condition at infinity |31[]. It can be shown that 5 is in fact the 
same function as that denoted C, in the fully relativistic treatment of axisymmetric stars of 
Ref. where the proper handling of this subtlety is described and an integral expression 
(Green's function expansion) for C, is given. The first term of this expansion can be adapted 
as the solution to Eq. ([17|) . 



III. HYDRODYNAMICS IN THE NEWTONPLUS APPROXIMATION 

The utility of the NewtonPlus approximation will be greatly enhanced if existing multi- 
dimensional hydrodynamics codes can be adapted to its use. As a concrete case we consider 
the Virginia Hydrodynamics code (VH-1) [Q, an implementation of the piecewise parabolic 



method This code treats a multidimensional problem with operator splitting: Succes- 
sive sweeps through the spatial dimensions are taken, with each sweep involving a Lagrangian 
step forward in time followed by a remapping of the fluid variables to a fixed Eulerian grid. 
By way of example, we here present the equations for a spherically symmetric calculation 
in the NewtonPlus metric, which are analogous to those solved by VH-1 in its Lagrangian 
time steps in a similar Newtonian calculation. Similar equations can be derived for use in 
the angular sweeps. 

The basic hydrodynamic equations can be expressed as conservation laws which follow 
from the vanishing divergences of the baryon flux vector and perfect fluid stress-energy 



tensor. It has been known for some time |^ that relativistic Eulerian hydrodynamics can be 
put into a conservative form amenable to methods developed for Newtonian hydrodynamics. 
In seeking to adapt VH-1 to an approximate relativistic treatment, we follow Appendix B 



of Ref. which deals with spherically symmetric relativistic hydrodynamics in Eulerian 
coordinates. It was shown that these Eulerian equations can be cast in a form that, while 
not truly "Lagrangian" in the sense of being in a comoving frame, is nevetheless similar to 
the Newtonian Lagrangian equations. 
Consider a vector of state variables, 

VL={D,S,E). (18) 
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The components of this vector are defined by 



D = TpB, (19) 
S^T\p + p)v, (20) 
E = T\p + p)~p, (21) 

where T = (1— f ^)^/^ is the boost factor between the fluid rest frame and the flxed "Eulerian" 
frame; ps, P, and p are respectively the baryon rest mass density, total energy density, and 
pressure in the fluid rest frame; and v is the fluid physical radial velocity as measured by 
an Eulerian observer. This state vector is governed by the hydrodjTiamics equations in the 
metric of Eq. (|), 

9iU + ^5,(rV-*f) = o-, (22) 



where the fluxes are 



and the sources are 



{Dv,Sv + p,S), (23) 



aD = 3Ddt<^, (24) 

"2 



as = -e^^^^{E + p) dr{^ + 5) + e^+^K 



p 



r 



+ dr{5 - $) 



+ 459*$, (25) 



(JE = -e'+"^S dr{<l> + 6) + {3p + Sv + 3E) dt<l>. (26) 

It is convenient to separate out the terms involving the advection of the state variables; to 
this end, Eq. (|2^) can be rewritten as 



dtu + — drir'e'-'^uv) + ^ dr{r^e'-''k) = a, (27) 



where the vector 

k={0,p,pv) (28) 

has been deflned. 

Next we consider how the conservation laws of Eq. ( P?]) can be treated by a hydrody- 
namics code like VH-1, which performs Lagrangian hydrodynamic time steps followed by a 
"remap" to an Eulerian grid. We begin with baryon number conservation. The baryon flux 
vector is = psu^, where is the fluid 4- velocity. The number of baryons in a proper 
three-volume described by the 1-form AS^ with edges Ar, AO, and A0 is 

constant = AS^, = (e'^^r^ sin 9 Ar A^ A^)TpB. (29) 

In VH-l, the mass density is updated during a Lagrangian step as follows. From the time 
averaged fluid velocities at the zone edges obtained from the solution of the Riemann prob- 
lem, the new positions of the zone edges are computed. From the new zone edge positions, 
the new zone volume is computed, and the new zone density is given by 
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(old zone volume) , . 

(new density) = (old density)- -. (30) 

(new zone volume) 

Clearly, Eq. (^) can be used in precisely the same way. There is a factor of F that must be 
dealt with, but presumably a relativistic Riemann solver |^ can be used if allowance for 
relativistic velocities is desired. 

In VH-1 the updates of the velocity v and specific internal energy e make use of the 
(Newtonian) Lagrangian equations 

dtV + r'^dmP = g, (31) 
dte + dmir'^pv) = v g. (32) 

Here the mass coordinate is defined by 

rr 

m = pbt"^ dr, (33) 



and 5^ is a force — gravity, for instance. 

In the NewtonPlus approximation a similar set of equations can be obtained. We trans- 
form from the variables t, r to t, m, where 

i=t, (34) 



m 



rTpBe-^''r^dr= De-^'^r^ dr. (35) 
Jo Jo 



Even though i = t, we still write derivatives as dt or dt to indicate whether m or r is being 
held constant respectively. We now use this change of coordinates to bring equation 
into a useful form. The first component of equation (p7|) is 



dtD + —dr{r^e^-^Dv) = 3D dt<^. (36) 



Subtracting u/D times Eq. (|36|) from Eq. ([27| ) yields 



Ddt{^'^+ De'^'%dr (^) + ^dr{r'e'-''k) = ct - 3ndt^. (37) 



Equations ( ^4[]3^) can be used to show that 

dt = dt-e'-''rhDd^, (38) 
a, = rV3*Da^, (39) 

so that Eq. { ^7\) can be expressed as 

d-t (^) + dUr'e'-^'k) = ^ - 3u9i$. (40) 
The second and third components of this equation are 

dt (^) + rV-" dmP = g + ^ (41) 
d, (I) + dUr'e'-^pv) = vg + 9,$, (42) 
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where the "gravitational force" g is defined by 



9 = CE + p) ^^^^ ^ ^ -e^-*(E + p) dU^ + 6). (43) 

In keeping with our practice of dropping exphcit time derivatives of $ (at constant r), the 
last two terms of Eqs. (^) and (^21) can be neglected. Furthermore, it is convenient^] to 
define a specific internal energy U hj U = E — D, so that on the left hand side one can 
write dt{E/D) = dt[{U/D) + 1] = diiU/D). 
The resulting equations, 

di (^) + r'e'-" dmP = 9, (44) 

dt(^^)+dm{r'e'-''pv)=v9, (45) 

are of the same form as Eqs. (|31|) and (|32|), allowing similar computational methods to be 
employed in their solution. It may seem surprising that the equations look so similar in 
the NewtonPlus and Newtonian formulations. As mentioned previously, it has been known 



for some time [^] that relativistic conservation laws can be cast in a "conservative" form 
similar to the Eulerian Newtonian equations. Because we have taken care to define a mass 
coordinate based on the proper relativistic 3- volume element, and chosen a metric with 
a close connection to the one giving the Newtonian limit, we have also been able to find 
equations quite close to the Newtonian Lagrangian formulation. 



IV. TESTING NEWTONPLUS GRAVITY WITH RAPIDLY ROTATING STARS 

In this section we present calculations of neutron stars undergoing rapid uniform rota- 
tion in order to assess the strengths and weaknesses of the NewtonPlus approximation to 
relativistic gravity. Our models were computed with a code described in Ref. which 
was written to compute the structure of relativistic axisymmetric stars. We have modified 
the code to include the ability to perform computations in the Newtonian and various New- 
tonPlus limits: with vanishing metric function 6, with "linearized" 6 (i.e. ignoring 6 on the 
right hand side of Eq. (^)), and "full" S (solving Eq. (|T^) as written). All of the Newton- 
Plus limits solve a two-dimensional (and stationary) version of the nonlinear Poisson-type 
Eq. d^) for the glorified "gravitational potential" $. For the results presented here, the 
high-density portion of the equation of state (EOS) is taken from Ref. [^, and is based 



on a field-theoretic description of cold dense matter. We also performed calculations with a 
polytropic EOS of adiabatic index 2, and found qualitatively similar results. 

Panel (a) of Fig. |l] shows mass vs. radius curves for spherical stars. The gravitational 
mass, a measure of total mass-energy, is defined for asymtotically fiat spacetimes by (see 
e.g. Ref. m 



Computationally, this prevents finite differences in the internal energy from being swamped by 
finite differences of the (rest+internal) energy. 
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M= {T\-T\)y/^dx^dx^dx^, 



(46) 



where here g is the determinant of the metric and are spatial coordinates. For the 
NewtonPlus metric, 

M = J e^-^^ {E + S\y sin 6 dr dO #. (47) 

From Eq. (^, it is evident that the proper equatorial radius is 

R = exp[-<l>(r = rsurface, = 7t/2)] r^, (48) 

where is the coordinate radius of the equatorial surface. Panel (a) shows that while the 
Newtonian limit exhibits no maximum mass with this EOS,[] the NewtonPlus approxima- 
tion does yield a maximum mass. Even with vanishing 6, the approximation captures this 
consequence of nonlinear gravity. The "linear 5" approximation follows the exact relativistic 
curve until the most dense configurations are reached, where the approximation "overshoots" 
the true mass (i.e. masses are too large by a significant factor when S is neglected, and a 
bit too small when the "linear 5" approximation is employed). Since pressure is the main 
source for 6 (see equation (|l^), the large pressures associated with such high densities raise 
6 to large enough values that it cannot be neglected on the right-hand side of equation (^7^ . 
As expected, the "full 6" approximation is indistinguishable from the relativistic results in 
spherical symmetry, where only two metric functions are needed to describe the spacetime 
exactly. The behavior of the "vanishing S" and "linear 5" approximations in comparison 
with the exact result is given further explanation in the Appendix. 

Panels (b)-(f) of Fig. |l| show various physical parameters of rapidly rotating configura- 
tions. In order to test the NewtonPlus approximation in a nonspherical setting, we ask the 
question: Given a definite number of baryons rotating at a given uniform angular velocity 
Q, what do the various treatments of gravity do with those baryons? (The fact that baryon 
number is a conserved quantity makes this an obvious way to compare different treatments 
of gravity.) To answer this question we have computed, for each treatment of gravity, a 
constant baryon mass sequence beginning at zero rotation (marked by squares) and ending 
at the mass shedding limit (marked by stars). The value of baryon mass chosen, 1.8 Mq, is 
close to the maximum baryon mass of 1.95 Mq for the equation of state we employed. The 
baryon mass is defined by (e.g., Ref. |]31|| ) 



^No turnover in the mass vs. radius curve appears in the Newtonian limit, up to the high-density 
boundary of the tabulated EOS. A configuration with central baryon mass density (total energy 
density) of 3.07 x 10^^ g cm^'^ (4.65 x 10^^ g cm~'^) has a gravitational mass of 15.6 Mq and 
radius 18.9 km in the Newtonian Hmit, while the relativistic configuration with this central density 
has a gravitational mass of 1.68 Mq and a radius of 9.40 km. We remind the reader that the 
Chandrasekhar mass phenomenon is a property of stars built on a polytropic equation of state 
with adiabatic index equal to 4/3, and that "realistic" nuclear equations of state will not generally 
exhibit this behavior. Instead, the upper mass limit of neutron stars derives from the the general 
relativistic instability indicated by the turning point in the mass vs. radius curve. 
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Mb = J {-n^Af')Vhdx^dx^dx^, (49) 
-^"^ TpBT^ sine drded(j), (50) 



where is defined after Eq. (^), is defined before Eq. (29), and h is tlie determinant 
of h^u = g^y + n^riy. Tlie second line (liere and in tlie equations below) is specialized 
to the NewtonPlus metric. The quantities plotted as a function of the (uniform) stellar 
angular velocity in panels (b)-(f ) are (b) gravitational mass; (c) equatorial radius; (d) angular 
momentum, given by (e.g., Ref. [pl[1 ) 

J = y T^^^/^ dx^dx'^dx'^, (51) 

= J e-^-^^ {E + p)Q sin^ e dr de d(j)] (52) 

(e) the eccentricity, defined here as 1 — rp/r^, where Vp and rg are respectively the polar and 



equatorial coordinate radii; and (f) the linear equatorial velocity, given by (e.g. Ref. ||3T[| ) 



t/=^(e^)X (53) 
= e-'-^'^rsin^ (54) 

where (e,^)^ is the basis vector corresponding to the coordinate 0. 

In panels (b)-(f), the efficacy of the NewtonPlus approximations can be judged by choos- 
ing a value of angular velocity and seeing how close the approximate quantities come to the 
fully relativistic value. While the "full S" approximation is indistinguishable from full rela- 
tivity in the spherical case, the two curves representing these treatments deviate from one 
another with increasing angular velocity. In panels (b) and (c), the "overshoot" of the "lin- 
ear S" approximation, discussed in connection with panel (a), is again visible at low angular 
velocities. At the highest angular velocities, the "linear 6" approximation actually appears 
to give slightly better results than the "full 5" calculations, but it must be remembered 
that this is a lingering result of the "overshoot" at low angular velocities. As expected, 
the angular velocities at mass shedding of the NewtonPlus approximations are closer to the 
relativistic values than the Newtonian case. The NewtonPlus treatments are quite successful 
at approximating the gravitational mass, radius, and eccentricity, while the success of the 
results for angular momentum and equatorial velocity is more modest. 

The relative success of the NewtonPlus approximation in determining these observables 
can be understood by examining expressions for these quantities in the fully relativistic case. 
Axisymmetric, stationary spacetimes can be described by four metric functions; the metric 
can be expressed pTl , |32 | 



g^pdx'^dx^ = -e^'de + e-2"GV2 sin^ e{d<P - N^dtf + ^^^-''\dr^ + rHe''). (55) 

In comparison, the NewtonPlus metric is conformally fiat (equivalent to requiring G = 
in Eq. (^)), and lacks the shift vector component N'^. In the relativistic case the angular 
momentum and equatorial velocity are given by 

J = J e^'^-^''G'\E + p){Q- N'^ysin^edrded(j), (56) 
U = e-^^Gr sin e{Q - N'''). (57) 
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In comparison with Eqs. ( p2D and (Q), an important difference is that is replaced by 
r2 — iV^. Since values of A^*^ can be a significant fraction of Q in extreme configurations (e.g., 
Refs. pl| , |39| ), the neglect of this shift vector component in the NewtonPlus approximations 
probably accounts for most of the difference in J and U with the relativistic case.0 The 
gravitational and baryon mass in the relativistic case are given by 

M = J e^^^-"^ G{E + S\ + 2e-''N^J^y sin 6 dr dO d(f), (58) 
Mb = J e^^-^"" G TpB sin 6 dr d9 dcj), (59) 

where = {E + p)e^'^Grsm6U, and the proper equatorial radius is 

R = (e-^GU,_,_,,=./2) r,. (60) 



In Eqs. (|59| ) and (§3), A^''^ does not appear at all, and it plays a minor role in Eq. 
This, together with the fact that the approximation of conformal flatness is known to work 



well for rapidly rotating stars shows why the NewtonPlus approximation gives good 
results for the gravitional mass, baryon mass, and equatorial radius. 



V. CONCLUSION 

Accurate neutrino transport, 3D hydrodynamics, and relativity are all essential for re- 
alistic supernova simulations. Given the constraints of current hardware, these cannot all 
be treated simultaneously with the detail they deserve. We have therefore presented an 
approximation to full relativity (or set of related approximations) that captures impor- 
tant relativistic effects in the quasispherical supernova environment: nonlinearity creating a 
deeper potential well, and pressure and internal energy density being nontrivial sources of 
gravitation. 

This "NewtonPlus" approach to gravity has a tight conceptual link with the Newtonian 
limit that yields certain advantageous features. The gravitational portion of multidimen- 
sional Newtonian calculations involves only the solution of the Poisson equation for the 
gravitational potential $, and taking its gradient to find the gravitational force. The basic 
idea of our NewtonPlus approach is to promote the Newtonian metric functions — which are 



^Since the lapse function and shift vector were quantities to be chosen at will in the (3+1) 
formalism as originally conceptualized, at first glance it may seem strange that the shift vector 
component N'^' should play a significant role in the determination of physical quantities. However, 
other coordinate choices have been made in arriving at Eq. (|55|), making the lapse and shift vector 
component quantities with physical content for which solutions must be found. Specifically, two of 
the four degrees of coordinate freedom have been used to choose the time coordinate t to label the 
spacelike hypersurfaces which are invariant under time translations, and the azimuthal coordinate 
<j) to measure the angle about the axis of symmetry. Once these two choices have been made, the 
metric components gtr, Qte^ 9^rj and g^e vanish (see Ref. |4C], which draws on ^^). The remaining 
two degrees of coordinate freedom have been used to set Qre = and gos = r'^grr- 
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linear in $ — to full exponentials. We also add a second metric function, 5, whose main 
source is pressure; hence this metric function vanishes in the Newtonian limit. The Ein- 
stein equations yield a nonlinear Poisson-type equation for a (now glorified) "gravitational 
potential," whose solution in 3D can be obtained in a manner similar to what would be 
done in the Newtonian limit. The inconsistencies in the Einstein equations arising from the 
reduced number of degrees of freedom turn out to be relegated to the subdominant met- 
ric function 6; they can be removed by ignoring angular variations in 6. This is expected 
to be successful in the supernova context because the region where 6 makes the greatest 
difference is where pressure is significant in comparison with rest mass density. Normally, 
this is the deepest portion of the collapsed core, which is roughly spherical even when the 
outer layers bulge at the equator due to rapid rotation.0 This strategy — allowing the main 
contribution to the gravitational field to be multidimensional and nonlinear, while allowing 
a spherical correction for the contribution of stresses — reproduces fairly well many of the 
physical characteristics of rapidly rotating (stationary) relativistic stars. In future dynamic 
calculations, we intend to ignore explicit time derivatives of $ that appear in the Einstein 
equations, and we have given specific arguments for the validity of this approximation in the 
supernova environment. (This can be tested in the future by comparing the results of spher- 
ically symmetric collapse simulations using NewtonPlus gravity with similar fully relativistic 
computations.) Importantly, the hydrodynamics equations in spherical symmetry obtained 
from the NewtonPlus metric are of the same form as those used in a popular Newtonian 
hydrodynamics algorithm, providing the expectation that existing Newtonian codes might 
be adapted fairly easily to the NewtonPlus approach. 

Finally, we make a few comments on approximations to full relativity employed by Math- 
ews and Wilson [45] and by Shibata, Baumgarte, and Shapiro |46] in binary neutron star 
merger calculations. While to our knowledge these approaches have not been implemented 
in supernova simulations, they conceivably could be; hence we offer a few comments by way 
of comparison. These approaches involve the use of a conformally flat metric, the imposition 
of a traceless extrinsic curvature tensor, and the neglect of explicit time derivatives of grav- 
itational variables. Some arguments in favor of these approximations are given in Ref. ||47||.p| 
These approximations are similar to ours in that both retain a reduced number of degrees 
of freedom (they keep four, while we keep two), both assume conformal flatness, and both 
ignore explicit time derivatives of gravitational variables. An important difference between 
Refs. and is that the latter neglect certain nonlinear gravitational terms that do not 



^Ultra-strong magnetic fields or differential rotation (e.g., Ref. |4^) can give rise to off-center 
density maxima, making the spherical correction for stresses via the metric function 6 less useful. 
The NewtonPlus approximation with 5 = could still be employed in such (probably exceptional) 
cases, however. 

^Mathews and Wilson's initial calculations produced the controversial result that the individual 
stars collapsed to black holes prior to merger, and it was suspected that this might be caused 
by their approximations rather than a physical effect. However, an error in their equations was 



pointed out in Ref. |48|, and the discrepancies with previous expectations and the results of other 



groups' calculations have been greatly mitigated in their most recent results [45| 
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arise in the limit of spherical symmetry. We would argue that an advantage of our approach 
is the conceptual link with the Newtonian limit. In the quasispherical supernova context this 
makes for only one dominant metric function — the glorified "gravitational potential" — while 
in the other approaches both the conformal factor and the lapse function are significant and 
must be solved for in multidimensions. They retain the entire shift vector; we have sug- 
gested dropping it altogether, though our work shows that this causes errors in the angular 
momentum and equatorial velocities. In the context of assumed conformal flatness, Wilson 
and Mathews also found that proper treatment of the shift vector was critical to accurate 
treatment of the angular momentum While angular momentum is perhaps less crit- 
ical in supernovae than in binary mergers, it is known that the assumption of conformal 
flatness alone (while retaining an azimuthal shift vector component) yields accurate results 
for rapidly rotating stars P2[. This suggests a future extension to the work presented here 



that would improve results for angular momentum and equatorial velocity. The azimuthal 
shift vector component iV^ obeys an elliptic equation, whose solution can be written as a 
Green's function expansion (see Ref. for an explicit formula). Like the solution of Eq. 
(p!7|) for 6, obtaining the leading term of A^"^ would involve only a pair of angular and radial 
integrations over the spacetime, but it would likely improve results for angular momentum 
and equatorial velocity significantly. 
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APPENDIX: 

This appendix shows the relationship between the NewtonPlus metric in the static limit 
and the (interior) Schwarzschild metric. We point out that the equation of hydrostatic 
equilibrium derived with the NewtonPlus metric is equivalent to the familiar Tolman- 
Oppenheimer-Volkov (TOV) equation, as expected if the NewtonPlus metric provides an 
exact solution in spherical symmetry. The connection to the Schwarzschild metric also al- 
lows greater insight into the behavior of the "vanishing 6" and "linear 6" approximations in 
the static limit. 

The Schwarzschild line element is given by 

ds^ = -^^de +(l-'^y\r^ + r^d^\ (Al) 



r 

where $ and m are functions of r. Comparing with a static version of Eq. (|^) in which $ 
and 5 are only functions of r, we find that equivalence between the two metrics requires 

r = e-*r, (A2) 

$ = $ + (5, (A3) 

'l+r— = 1+r— . (A4) 
dr j \ dr j 
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We now describe how equations derived from the NewtonPlus metric are equivalent to 
the TOV equation, 



dr r \ r J \ r 

The rr component of the Einstein equations in the static NewtonPlus metric is 

2 



(A5) 



2d6 ^d6d<^ 
r dr dr dr \ dr . 



ine p. 



(A6) 



Using Eqs. 
equation becomes 



I), and 



it is a straightforward exercise to show that the TOV 



dp 
dr 



, (d^ d6^ 
\dr dr . 



(AT) 



which is the static version of Eq. (|1|) (see also Eqs. (|3|) and (H)). 

It is well known that m(r) evaluated at the stellar surface is equal to the gravitational 
mass. Using Eqs. (|A4|) and (|A6|) , we find 



m 
r 



d^ 



dr 
dr 




2d6 ^d6d<^' 
r dr dr dr , 



^ -2$ 

ine p 



(A8) 
(A9) 



By evaluating these expressions at the stellar surface, we gain insight into the behavior of 
the "vanishing 6" and "linear 6" limits. In the exact solution, the second term in Eqs. 
( [Asp and ( |A9| ) is negative. In the "vanishing 5" approximation, this term becomes positive, 
explaining why the gravitational masses obtained in this limit are too large. Because the 
exact second term is negative, the quantity in parentheses in Eq. ( [A9| ) must be positive and 
greater in magnitude than the pressure term. Also, since 6 is subdominant, Eq. ( [A^ ) shows 
that d6/dr is positive. The "linear 6" approximation involves neglecting terms containing 
6 that are nonlinear in metric functions, which here involves neglecting 2{dS/dr){d^/dr). 
This makes the quantity in parentheses in Eq. (^^) more positive, which makes the overall 
second term too negative. This explains the "overshoot" (gravitational masses too small) 
seen this limit. 
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FIGURES 



FIG. 1. Panel (a): Mass vs. radius curves for spherical configurations computed with various 
treatments of gravity. Panels (b)-(f): Various physical quantities characterizing uniformly rotating 
configurations, plotted as a function of angular velocity. Each curve represents a constant baryon 
mass sequence computed with the treatments of gravity labeled in panel (a). 
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